function VAR = proxy_svar(VAR,bias)

if nargin==1
    bias = 0;
end

% Preliminaries
data_Y = VAR.data.Y;
data_Z = VAR.data.Z;
 
[T, n_y] = size(data_Y);
[T_z, k] = size(data_Z); 
 
p = VAR.p; 
T_y  = T - p;  
[overlap,overlap_idx] = intersect(VAR.data.time_y,VAR.data.time_z,'stable'); 
ty_zstart = overlap_idx(1);
ty_zend = overlap_idx(end);

% Estimate
X          = lagmatrix(data_Y,1:p);
X          = X(p+1:end,:);
Y          = data_Y(p+1:end,:);
VAR_coeff  = [X ones(length(X),1)]\Y;
VAR_res    = Y-[X ones(length(X),1)]*VAR_coeff;
VAR_res_restr = VAR_res(ty_zstart:ty_zend-p,:);
Sigma      = (VAR_res'*VAR_res)/(T_y-n_y*p-1);
Sigma_restr  = (VAR_res_restr'*VAR_res_restr)/(ty_zend-p - ty_zstart+1  -n_y*p-1);

% collect results
VAR.T      = T_y;
VAR.p      = p;
VAR.k      = k;
VAR.n      = n_y;

% VAR.beta   = VAR_coeff; 
VAR.beta   = VAR_coeff + bias; 

if nargin>1
     % stability check
    bb = zeros(VAR.n,VAR.n*VAR.p);
    bb(1:VAR.n,1:VAR.n) = VAR.beta(1:VAR.n,1:VAR.n);
    for jjj=2:VAR.p
        bb(1:VAR.n,1+(jjj-1)*VAR.n:jjj*VAR.n) = VAR.beta(1+(jjj-1)*VAR.n:jjj*VAR.n,1:VAR.n);
    end
    d = stability(bb,VAR.p); 
    if d >= 1; error('--> WARNING: the estimated VAR model is not stable!'); end
end


VAR.Sigma  = Sigma;
VAR.Sigma_restr  = Sigma_restr;
VAR.res    = VAR_res;
VAR.m      = data_Z(p+1:end);
VAR.ty_zstart = ty_zstart;
VAR.ty_zend   = ty_zend;

% Identification: only on restricted sample  
Phib = [ones(length(VAR.m),1) VAR.m]\VAR_res_restr;
Phib = Phib(2:end,:);
Phib11  = Phib(1:VAR.k,1:VAR.k);
Phib21  = Phib(1:VAR.k,VAR.k+1:VAR.n);
b21ib11 = (Phib11\Phib21)';
Sig11   = VAR.Sigma_restr(1:VAR.k,1:VAR.k);
Sig21   = VAR.Sigma_restr(VAR.k+1:VAR.n,1:VAR.k);
Sig22   = VAR.Sigma_restr(VAR.k+1:VAR.n,VAR.k+1:VAR.n);
ZZp     = b21ib11*Sig11*b21ib11'-(Sig21*b21ib11'+b21ib11*Sig21')+Sig22;
b12b12p = (Sig21- b21ib11*Sig11)'*(ZZp\(Sig21- b21ib11*Sig11));
b11b11p = Sig11-b12b12p;
b11 = sqrt(b11b11p);
VAR.b1 = [b11; b21ib11*b11];
% VAR.b1 = VAR.b1(:,1)/VAR.b1(1,1);  % apparently not normalized by GK

% Structural shocks
% VAR.b1 = [ 0.1955
%    -0.0328
%     0.0289
%     0.1130];
VAR.eps1 =  VAR.res * VAR.b1;
 
% Impulse Responses 
irs(VAR.p+1,:) = VAR.b1;
for jj=2:VAR.IRF_hor
    lvars = (irs(VAR.p+jj-1:-1:jj,:))';
    irs(VAR.p+jj,:) = lvars(:)'*VAR.beta(1:VAR.p*VAR.n,:);     
end
VAR.irs = irs(VAR.p+1:end,:); 
 
end